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specific heat at constant volume and constant 
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internal energy, Eqn (3) 
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Single crystals are solid with the most uniform atomic 
structure that can be attained. This uniformity is the basis for 
a wide range of technological applications of single crystals. 
The most important feature of crystal is that it allows the 
transmission of acoustic, electromagnetic and charged particle 
with a minimal scattering. In particular, this transparency can 
be selectively modified by the controlled addition of 
impurities (dopants) to produce semi-conductor materials. [ 1 ] 

Crystals are made from carefully controlled phase changes 
that involve solidification of a high temperature melt (melt 
growth), growth from low temperature mixture (solution growth) 
and condensation of a vapor (vapor growth). In all of these 
methods, the the nutrient materials are in either liquid or vapor 
phase. The physical transport mechanism at the interface between 
a nutrient fluid and the solid surface of growing crystal is 
controlled by mass diffusion. The solute gradient at the 
interface is , however, critically dependent upon the solute and 
thermal distribution of the nutrient fluid far away from the 
interface. This far-field flow behavior is controlled by the 
convection which is produced by various physical mechanism. 
Convection can be generated by the density gradients in a 
gravitational field (natural convection ), or by an externally 
applied pressure (forced convection) or by both (mixed 
convection). In a low gravity environment, surface tension driven 
convection can be important when the nutrient fluid is allowed to 
have a free surface. Thus, in general, the growth of crystals is a 
coupled process of heat and mass transport, fluid flow , phase 
transformation and chemical reactions. 

Even though the technology of crystal growth is well grounded 
in physical chemistry and materials science, the treatment of 
transport processes sometimes have been rudimentary and 
qualitative. To improve the quality of crystals, it is apparent 
that the flow characteristics at the growing interface must be 
understood so that one can control whatever influences the 
convection might have on the quality of crystals. Significant 
attempts have been made in recent years by fluid dynamicists to 
understand the complex interactions among various transport and 
physical processes in crystal growth. [1,2, 3, 4] 


II. OBJECTIVES OF THE PRESENT RESEARCH 

Much of the earlier works on the study of convection in 
crystal growth were based on the premise that convection is 
deleterious and therefore should be kept minimum at all cost. 
Materials processing in a space environment would provide a low 
gravity environment which will reduce buoyancy induced natural 
convection significantly. It is not likely however to create a 
completely acceleration free environment. A space station in 
which materials processing would occur is subject to microgravity 
due to near earth orbit and viscous drag force. Diffusion 
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in crystal growth are usually extreme— y s_ow an- 
natural convection created by transient microgravity _ could be in 
the same order of magnitude. Surface tension driven natuiax 
convection, which is' usually insignificant in terrestrial 
conditions, becomes very important in low gravity environment 
when the crystals are made by a containerless process. '23 


. he 


purpose of the present research is to initiate a 
development of a numerical model which can be used to simu-ate 
crystal growth in realistic space environments. The thermal and 
solutal fields encountered in crystal growtn from; solution are, 
in general, three-dimensional and time-dependent. Transport- 
properties, such as thermal diffusivity, mass diffusivity and 
fluid viscosity are i general functions of solute ^concentration 
and fluid temperature of the mixture. Solidification ^ processes 
usually involve liberation of latent heat at the interface. Also, 
rate of solidification depends on segregation 
coefficient , which in turn depends on f * u i u .emperat 
scale three-dimensional moce- '•■•’hich accrur.os^ - - - 
physical conditions is not warranted until sufficient 
of simplified two-dimensional situation 


* str i but ion) 
ure . A full 

1 “t h 0 S 0 

analyses 


ire performed. 


a first ■=. ten , a two-dimensional, time— dependent naturaj 
convection driven* by both thermal and solutal buoyancy forces is 
selected in this research . Sven though this physical mcdex is 
far simpler than the actual processes, it should provide insight: 
into the complex interactions between thermal and soiUvi. 
convections encountered in crystal growth. Detailed description; 
on the physical model and appropriate equations are presentee i: 
section III. Numerical results are presented in section 
followed by conclusions and recommendations m section 
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Two-dimensional thermal and solutal convection ^ across^ a 
rectangular cavity is uepxu^eu^ 
are maintaint-ed at isothermal 
component A and the horizontal. 
and impermeable to mass diftusi 
made of two components of liqui 

molecular weights. Heat transfer occurs from the 

temperature surface (left wall > to the -ow temperature s-r- a- 
(right wall). Mass transfer of species A also occurs from 
’•eft surface where the mass concentration of A is higher tha 
the right surface. It is assumed that all surfaces 

impermeable to species B. One might consider species a tc^e 
dopant in melt growth or crystal material m vapor or solution 
growth. Relative location of higher temperature wa-- an: -g.^r 
concentration surface is immaterial. 
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Thermal conduction occurs from the hot surface to the r - al y 
and the temperature of the fluid increases over its mitia- 
value. Since the volume 


most fluid expands as increasing 
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temperature, density of the fluid becomes less than its 
surrounding fluid. Lighter fluid elements rise up when the 
gravity acts along downward. This buoyancy force sets the 
convective flow motion within the cavity. Buoyancy force will act 
along the downward direction next to the cold surface at the 
right. Unless the initial temperature difference between the 
walls and fluid is large, one amount of heat transfer is small 
and the density change is almost negligible. This is known as the 
Boussinesq approximation. Mass diffusion from the left surface 
to the fluid will have similar effects on the flow field. Suppose 
species A has higher molecular weight than species B. Then 
diffusiuon of species A into the fluid will increase the density 
of the mixture and the soiutal buoyancy force will act downward 
near the left wall. If the molecular weight of species A is lower 
than E, the soiutal buoyancy force will be in the opposite 
direction. Since the mass diffusion through liquid is extremely 
slow, Boussinesq approximation can be invoked for mass transfer 
as well. 


Governing Equations 

The governing equations for the conservation of mass, 
momemtum, energy and species A can be expressed by a generalized 
transport equation of the form 



Table 1 contains appropriate expressions for (j> , p and S for 
each equations. 

Buoyancy forces due to thermal and soiutal convections are 
included in the vertical direction momentum equation. All 
transport coefficients are assumed to be constants. No latent 
heat liberation is considered. The pressure appearing in the 
momentum equation represents only dynamic pressure. Static 
pressure is eliminated by using the density variation given by 
the Boussinesq approximation, i. e.. 


/>Y°[ 1 ~ ^ ^ ~ w a > ~ f 3 (T-~r°) J. 




It is assumed that the mass concentration of species A of the 
liquid at the right vertical wall is equal to the mass 
concentration of species A at the solid. This assumption amounts 
to having a unit segregation coefficient at the growing surface. 
It is also assumed that the interfacial mass flux of species A at 
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C £ 1 C U _i_ 3. t • 0 Cl 


:-ne mas* 


averaged 


.ow 
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in the horizontal direction should be imposed at these walls for 
a more accurate boundary treatment.) Internal energy of the fluid 
is aiven bv 
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Initial and Boundary Conditions 

Initial and boundary conditions depend on the physical 
problems under consideration. Two types of initial conditions are 
considered; one corresponds to a static equilibrium condition 
and the other corresponds to a steady-state thermally driven 
natural convection. For the static initial conditions 

u y, *, 3, -t= o j — o 

/° V ) * p* 

T t-o = T ‘ <4) 

W A = VU4 0 


Boundary conditions are fixed for all numerical results 
reported in this paper. They are 


T 3, t) - T» , Tfx=*v », t) = 

* W <(l / hfc * M0ie 


1 

1 _ J 



= » / 

X, V-o, t 



Nondimensional Parameters 


The present numerical method utilizes primitive variables in 
dimensional forms. It is convenient f however, to introduce 
appropriate flow characteristics into the governing equations to 
obtain important nondimensional parameters which govern the flow. 
Numerical and experimental data presented in nondimensional 
parameters are universal and can be used to infer other similar 
problems with different fluid properties and physical geometries. 

Nondimensional variables are defined by 
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U. 




T - 


t-t° 


AT 



A\* 1 a 


(6) 


introducing these nondimensional variables into the governing 
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equations (1), it can be readily shown that nondimensional 
equations contain five nondimensional parameters. They are; 
aspect ratio (a), thermal Grashoff number < Gr ), solutal Grashoff 
number (Gr ), Prandtl number (Pr), and Schmidt number (Sc), 
Numerical results will be presented in terms of nondimensional 
quantities . 


VI. ANALYSIS 


Numerical Method 


The numerical method used for the present study is 
volume method based on the SIMPLE algorithm. [5] Since a 
description on the essence of SIMPLE can be found in a 
[5], it will not be discussed in this report. 


a finite 
detailed 
ref erer.ce 


The SIMPLE method, in principle, can be applied to both high 
speed flow as well as low speed flow by using a proper equation 
of state and including the density variations due to the pressure 
change. These modifications were added to the original SIMPLE 
method and the revised method were successfully applied to a 
transient acoustic wave propagation in a water hammer 16] and a 
transient natural convection in a square cavity [7], 


Physical and Numerical Parameters 

Physical and numerical parameters used for the present study 
are listed in Table 2. Physical properties are selected only to 
give representi ti ve thermal and solutal Grashof numbers, Prandtl 
and Schmidt numbers encountered in solution growth under 
terrestrial conditions. Thus these values may not represent any 
actual crystal growth environment. Numerical parameters are 
selected (by trial and error) such that either the solution 
approaches to the true solution (if exists) or converges to 
steady or quasi-steady solution within a reasonable computing 
time . 

Numerical Results 


Four different cases 
effects of thermal and 
interactions . 


are simulated to investigate the 
solutal convections and their 


CASE Is Thermal convection : Nondimensional parameters used for 

this case are 


p 

£r T = = O ; pr = 7®7 ) 4 = o.r 3 

£ x io c 

This case is studied first to check the accuracy of the 
numerical code developed. Thermally driven natural convection in 
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rectangular geometries are extensively studied and the present 
results can be easily checked against existing data. 

Numerical accuracy and computational economy depend on the 
choice of computational mesh arrangement. Since the driving force 
for the present convection comes from the vertical walls. The 
resulting buoyancy driven convection will form thermal and 
momentum boundary layers. It is therefore necessary to have fine 
grids near the vertical walls. Invoking scale analysis for the 
boundary layer regime [8], the thermal boundary layer thickness 
is 

rs ~'A - 3 

d T ~ H R<K T *• A. *7 K 16 ** 1 . 

11 ( 8 ) 

The .-.ash size at the vertical walls therefore should be small 
enough that at least- one computational mesh should be placed 
within this boundary layer. The mesh arragement along the x- 
di recti on is given by 

<5X„ r 0+ -“Vr > (9) 

r where <TX| is “the first mesh size and AXftr is the percent 

increment. The mesh arrangement near the horizontal walls need 
not be as small walls since the horizontal walls are adiabatic 
and the momentum boundary layer is thicker than the thermal 
boundary layer ( Pr > 1.0), Grid expansion along the y-direction 
is similar to x-direction. Number of trial runs were made before 
settling to the mesh arrangement as shown in Figure 4-a- Total 26 
(X-direction) by 20 (Y-direction) meshes were used in all 
subsequent computations. (Actually, the first mesh size near the 
vertical walls are much larger than the thermal boundary layer 
thickness. However, solut-al boundary layer thickness, discussed 
later, dictates this size.) 


Figure 2 shows transient behavior of rate of energy transfer 
across the cavity evaluated at the center of the cavity. Nusselt 
number is defined by 


f\i W £ = 




H AT 

L 


( 10 ) 


where (SJ represents the total energy flux given by 

n 


Q 


H* f [/ ,uc r( T - T °; ^ 

a ' y = i 


(11 ) 
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Nusselt number increases rapidly as the hot and cold fluid 
intrude the top and bottom portions of the cavity • Energy 
^^®nsfer is at the maximum when the convection velocity reaches 
the maximum at this early stage of transient convection. 
Numerical solution exhibits oscillatory horizontal flow speed. 
These oscillations are mainlyresponsible the oscillatory Nusselt 
number. Oscillatory energy transfer across a square cavity was 
theoretically predicted by Patterson and Imberger [9] and was 
numei ically confirmed. [7,9] Oscillation periods «n d amplitudes 
for the present case (a=0.53> are different from those for a=l . 0 
reported in ref [7, 9] . 

Steady- state Nusselt number computed by the present 
simulation is about 27. Bejan [8] obtained the analytical 
solutions for the boundary layer regime with Pr>l. Figure 3 
shows the average Nusselt number obtained by his solutions. Using 
his results , the average Nusselt number for the present case is 
about 29, which is in reasonable agreement with the computed 
value . 



0.364 


Figure 3. Nusselt number correlations [8] 


Figure 4 shows the global view of steady-state flow and 
temperature fields for Case 1. Very thin thermal and momentum 
boundary layers are developed at the vertical walls. Temperature 
in the core regions are highly stratified and isotherms are 
almost parallel. Most of the fluid in the core is relatively 
stagnant. 

Figure 5 shows the temperature and velocity profiles at two 
selected spatial locations in the cavity. Figure 5a shows the 
vertical velocity and fluid temperature at a horizontal plane 
located at the mid— height of the cavity. Thermal boundary layer 
is located within the momentum boundary layer and the over-all 
boundary layer thickness is in the order of % of the cavity 

width. Figure 5b shows the horizontal velocity and temperature of 
the fluid at the vertical plane located at the mid-width of the 
cavity. Magnitude of horizontal velocity and temperature 
*ii s ^ r ibution at this plane as shown in Figure 5b determines the 
energy transfer rate across the cavity. Note that the vertical 
velocity is order of magnitude higher than the horizontal velocity. 
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CASE 2; Solutai convection ; Nondimensional parameters used for 
this case are 

- J,3 xio ff > <^<- T - o ; S* C = 2 I 40 ; fl=o,r3 

» L7 * XI o to * - 2 * 

Nondimensional solute concentration equation is identical to 
the nondimensional energy equation except the Schmidt number 
appears in place of the Prandtl number. By virtue of this 

analogy between equations, average mass transfer of species A 
across the cavity can be determined by using the Nusseit number 
correlations with simple replacement of nondimensional parameters 
such as 


Pr 4 Sc 

Gr T 

* 

Nu * 

Ck 


*13) 

the results 

of 

Nusseit 

number 

correlation 

a= 

shown in 


Figure 3, the average steady-state Sherwood number corresponding 
to the given mass Rayleigh number is approximately 226. 


Figure 6 shows the computed transient Sherwood number 
evaluated at the center of the cavity defined by 


Q 


M 


= 


where is the net mass transfer of species A at the center of 

the cavity and given by 

< 2 „-- J -fhs — (15i 

Mass transfer increases sharply as the convection current reaches 
the center of the cavity and undergoes a weak oscillation before 
approaching the steady value of about 222. This is m exceiien- 
agreement with the analytical result. 

Using the analogy, solutai boundary layer thickness can be 
computed as 


( ± H > 


£ 


_ -'A -<f 

H 'X' 2 , n XI O 


This implies that 
vertical wall shou 
system =0.0005) 
size of the control 
( <Ty, =0.001) the 
Sh=97 . 


the center of the control volume next to the 
Id be less than 2.S4E-4 m. The present grid 
does satisfy this requirement. When the grid 
volume next to the vertical walls was doublea 
steady Sh number was 17Q and when OX | = 0.vu2, 
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The amplitude and frequency of Sherwood number oscillations 
are much smaller than the Kusseit number oscillations as shown in 
Figure 2. This is due to extremely high Schmidt number (2140' 
used in Case 2 compared with a smaller Prandtl number (7.0) 
used in Case 1. Viscous damping is more pronounced in Case 2. 
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"ions have double humps. Mass transfer across the 
is however mainly contributed by the convection current 
the horizontal walls since the mass concentration remains 
at the initial ■ . . ne around the second velocity maximum. 


Layered convec t i o n 
water < saturated with 
transfer to water ccc 


could occur, for example, when a stable 
salt) is heated from the side wall. Heat 
rs from the wall by conduction . Heated 
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Figure 7. Steady-state velocity vectors (a), stream lines (b) and 
iso-concentrations for Case 2 



water expands and rises up to the level where its density matches 
the surrounding density and flows horizontally away from tne net- 
wall. If the mass diffusivity of salt through water is equai to 
the thermal diffusivity, density chance cue on heat transfer -u_ 
mass transfer can cancell each other. -.n this case.^ _ a y e r e d 
convection can not occur. Double— dir fusive mow m vcivt= t-wu 
distinct diffusion processes with different diffusion 
coefficients. The present results are somewhat similar to go.id-.g- 
diffusive flow situation even though mass transfer is the only 
available diffusion mechanism. No satisfactory ^ physical 
explanation can be offered at present moment to explain tne 
layered convection in Case 2. 


CASE 3. Opposing combined convection : — this case both thermal 
and solutai buoyancy effects are considered. The mitia^ f« ow an ^ 
temperature fields are those obtained m _ase 3 ° m m ^ buoyancy 


force acts 


he opposite cire^ 


t h e the rm a 1 


'orce. Nondimensional parameters are 


<v- T =q*/o r ; 3x/ ; ?r - 17 ' 07 ^ 

a~o,r3) R* T = 6\36 Xio* ; ^ = i.7fxio'° <17; ' 

Two additional nondimensional parameters are defined in order _ to 
compare the thermal and solutai effects. They are _e wis .i-sr.^er 
and the buoyancy ratio 
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Figure 9. Transient Nusselt and Sherwood 
number variations for Case 3 
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number for Case 1 and Sherwood 


steady state Mussel t number for case x ana snerwooa number .or 
Case *2 are also indicated for comparison. Due to the opposing 
buoyancy effects, energy and mass transfer rates are markedly 
dropped compared with a pure thermal and a pure sclutau cases. 
Numerical solutions at the end of computations were not steady- 
state. Small fluctuations in Sherwood number are due _ to 
fluctuating velocity field. Qualitatively similar fluctuations 
were observed in experiments 


r i i 


Figure 10 shows "the quasi “•steady state f 


concentrat ion x ^ _✓ •-* •- - * ■ »- — * — — 

loops can be identified: a slow counter clockwise circulation near 
the solid walls due to solutal effects and a dominant clockwise 
circulation in the rest of the cavity. Overall flow 
characteristics are determined by the thermal effects because of 
high Lewis number. Temperature inversion (hotter fluid on top of 
cold > occurs at the lower left and top right corner of the cavity 


at the end of computation 


:>w , temperature and 
Two convection 


10c> due to solu‘ 


0 f f 0 c t' = . 


(Fig 

CASE 4. Aiding combined convection : In this case, soxutal 

buoyancy force and the thermal buoyancy force are acting along 
bhe same direction. Nondimensional parameters are identical to 
those of Case 3 except the change i n sign 
Initial conditions are again provided by 


for the buoyancy ratio, 
ase 1 . 


Figure 11 shows the transient heat and mass transfer rate 
across = the cavity for the aiding case in comparison with a pure 
thermal and a pure solutal cases. Numerical solutions , _ in 
particular the mass transfer rate, exhibit the mild fluctuations 
toward the end of computations as Case 3. Quasi-steady-state 

Lding case increased markedly over Case 
number is decreased below that of Case 1. 
t is a consequence of thermal and solutal 
expained in a later section. 


Sherwood number for the 

2. However, 

Mussel t 

This unexpec 

ted resul 

interactions , 

which is 

Figure 12 

shows “the 


concentration fields for 


ruasi 
?ase 3 


-steady velocity, temperature and 
Over-all features are similar to 
and Case 2 (Figure 12d). 


case 1 ( Fiures 12a, 12b and 12c) 

Thermal and Solutal Interactions 

y q anal vse the therm ial and soiutsi interactions more closely , 
detailed flow profiles at a horizontal plane located at the mid- 
height of the* cavity are presented in Figure 13. It can be seen 
that opposing solutal buoyancy forces result in the negative 
vertical velocity near the left vertical wall and reduced 
vertical velocity compared with the pure thermal case. 
Temperature and concentration gradients at the left wall are also 
reduced compared with pure convection cases.- 

, vertical velocity is increased near the 
the solutal buoyancy force as expected. 


For 


t-.hi 


:-he aiding case 

vertical w a -l. -l due to ~ — j * - 

However, vertical velocity is actually smaller than that of Case 
1 in most of the remaining velocity boundary layer. Temperature 
gradient at the wail remains almost identical to Case 1. But 
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Figure 10. Quas i -steady-state velocity vectors (a), stream lines (b) , isotherms 
and iso-mass concentrations for Case 3 
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concetraticn gradient- is much higher t-han t-hat of Case 2 due to 
the increased vertical flow speed compared with Case 2. 
Consequently, Nusselt- number for this aiding case is slightly 
lower than Case 1 but the Sherwood number is much higher than 
Case 2. 

The velocity profiles shown in Figure 13a are similar to 
similarity solutions [11] over a vertical flat, plate, except the 
aiding case. (Figure 14) Direct one to one comparison is not 
possible since two problems differ in geometry and nondimensional 
parameters . 


(Sk t ) 
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isotherms (c) and iso-mass concentrations (d) for Case 
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Figure 13. Vertical velocities (a), temperatures (b) and mass concentrations (c) 
distributions in boundary layers at the left vertical wall at the mid 
of the cavity 
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Figure 15. Horizontal velocity (a), temperature (b) and mass concentration distributions 
at a vertical plane located at the mid-width of the cavity 
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Figure 15 shows the horizontal velocity, temperature and 
concentration profiles at a vertical plane located at the mid~ 
width of the cavity. These results again show that the mass and 
heat transfer across the cavity are strongly dependent on the 
interactions between the thermal and solutal convections and that 
opposing and aiding effects can have unexpected influences on the 
transport mechanisms. 

Effects on Crystal Growth 

Figure IS shows the temperature and concentration at the 
center of the control volumes next to the right vertical wall, 
which might represent the growing surface of a crystal. Highly 
nonuniform temperature and concentration profiles there imply an 
uneven heat transfer and crystal growth rate. As expected, 
convection enhances the rate of crystal growth significantly over 
the pure diffusion cases but with detremental effects on the 
crystal structure. 


T *-0,5 — 0. 3 S' | W/*-“0,F 



Figure 16 . Temperature and mass concentration profiles along the 
vertical coordinate at the control volumes next to the 
right vertical surface 


XIX-23 


V. CONCLUSIONS AND RECOMMENDATIONS 


Numerical 
convections in 


model based on the 
chosen to represent a 


simulations of thermal and solutal buoyancy^driven 
a rectangular cavity are 


SIMPLE algorithm, 
possible crystal 


performed by a numenca. 
Physical parameters are 
growth from solution. 


Numerical results for a pure thermal convection were found to 
agree reasonably well with a known solution. turner icax 
for a pure solutal convection also agree well wiv*« 
solution in terms of the average Sherwood number. Ho^er, 
velocity field for the pure solutax convector. , e '^ ibl “ s 
peculier feature which can not be explained. Relatively xarge 
solutal Reyleigh number (1.G7E+10) and an extremely 

number (2140) combination might have produced the ^np^/s^a. 
solutions. This has to be checked more carefully c.roug.. — 
numerical ex per i mentations . 

Overall energy and mass transfer rates are seen to depend 
strongly on the interactions between the^ thermal an-. s^-ai 
effects. These interactions are high-y nonlinear and t.. 
state conditions may not exist for certain comomea convec-iwu-. 
To arrive at the more conclusive and quantitative Correia , 

however, exes;? ioIls covering * 

nondimensional parameters are required. 

To investigate the double-diffusive effects on the flow, * 
non-uniform initial solutal distribution in the cavity^ is nee-J- 
before the commencement of heat transfer from siue waxa.s. = 
initial solute concentration was assumed constant for pie=ent 
analysis. ) This can be easily accomplished by assigning any non 
uniform solute distribution in the cavity. 


Af-ker performing the recommended 

verifications as discussed above, one can 
on the problem formulations, such as the i 
segregation effects and the temperature boun 


refinements and 
relax the constraints 
nterf acial velocity, 
dary conditions. 
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